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Quantum cavities or dots have markedly different properties depending on whether their classical 
counterparts are chaotic or not. Connecting a superconductor to such a cavity leads to notable 
proximity effects, particularly the appearance, predicted by random matrix theory, of a hard gap in 
the excitation spectrum of quantum chaotic systems. Andreev billiards are interesting examples of 
such structures built with superconductors connected to a ballistic normal metal billiard since each 
time an electron hits the superconducting part it is retroreflected as a hole (and vice- versa). Using 
a semiclassical framework for systems with chaotic dynamics, we show how this reflection, along 
with the interference due to subtle correlations between the classical paths of electrons and holes 
inside the system, is ultimately responsible for the gap formation. The treatment can be extended 
to include the effects of a symmetry breaking magnetic field in the normal part of the billiard or an 
Andreev billiard connected to two phase shifted superconductors. Therefore we are able to see how 
these effects can remold and eventually suppress the gap. Furthermore, the semiclassical framework 
is able to cover the effect of a finite Ehrenfest time, which also causes the gap to shrink. However 
for intermediate values this leads to the appearance of a second hard gap - a clear signature of the 
Ehrenfest time. 

PACS numbers: 74.40.-n,03.65.Sq,05.45.Mt, 74.45. +c 
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I. INTRODUCTION 

The physics of normal metals (N) in contact with su- 
perconductors (S) has been studied extensively for al- 
most fifty years, and in the past two decades there has 
been somewhat of a resurgence of interest in this field. 
This has mainly been sparked by the realization of ex- 
periments that can directly probe the region close to the 
normal-superconducting (NS) interface at temperatures 
far below the transition temperature of the superconduc- 
tor. Such experiments have been possible thanks to mi- 
crolithographic techniques that permit the building of 
heterostructures on a mesoscopic scale combined with 
transport measurements in the sub-Kelvin regime. Such 
hybrid structures exhibit various new phenomena, mainly 
due to the fact that physical properties of both the super- 
conductor and the mesoscopic normal metal are strongly 
influenced by quantum coherence effects. 

The simplest physical picture of this system is that the 
superconductor tends to export some of its anomalous 
properties across the interface over a temperature depen- 
dent length scale that can be of the order of a micrometer 
at low temperatures. This is the so-called proximity ef- 
fect, which has been the focus on numerous surveys; both 
experimental^— and theoreticaliSr— . 

The key concept to understand this effect^— is An- 
dreev reflection. During this process, when an electron 
from the vicinity of the Fermi energy (Ep) surface of the 
normal conductor hits the superconductor, the bulk en- 
ergy gap A of the superconductor prevents the negative 
charge from entering, unless a Cooper pair is formed in 
the superconductor. Since a Cooper pair is composed 



of two electrons, an extra electron has to be taken from 
the Fermi sea, thus creating a hole in the conduction 
band of the normal metal. Physically and classically 
speaking, an Andreev reflection therefore corresponds to 
a retroflection of the particle, where Andreev reflected 
electrons (or holes) retrace their trajectories as holes (or 
electrons). The effect of Andreev reflection on the trans- 
port properties of open NS structures is an interesting 
and fruitful area (see Refs. [l7lri8l and references therein 
for example), though in this paper we focus instead on 
closed structures. Naturally this choice has the conse- 
quence of leaving aside some exciting recent results such 
as, for example, the statistical properties of the conduc- 
tance 19 , the magneto-conductance in Andreev quantum 
dots 2 ^, resonant tunneling^! and the thermoelectrical ef- 
fec t 2 ^ 23 in Andreev interferometers. 

In closed systems, one of the most noticeable mani- 
festations of the proximity effect is the suppression of 
the density of states (DoS) of the normal metal just 
above the Fermi energy. Although most of the exper- 
imental investigations have been carried out on disor- 
dered systemsi^^, with recent technical advances in- 
terest has moved to structures with clean ballistic dynam- 
ics 2 -^ ^ 24 ' 25 . This shift gives access to the experimental 
investigation of the so-called Andreev billiard. While this 
term was originally coined 26 for an impurity-free normal 
conducting region entirely confined by a superconducting 
boundary, it also refers to a ballistic normal area (i.e. a 
quantum dot) with a boundary that is only partly con- 
nected to a superconductor. The considerable theoretical 
attention raised by such a hybrid structure in the past 
decade is related to the interesting peculiarity that by 
looking at the DoS of an Andreev billiard we can deter- 
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mine the nature of the underlying dynamics of its classi- 
cal counterpart 27 . Indeed, while the DoS vanishes with 
a power law in energy for the integrable case, the spec- 
trum of a chaotic billiard is expected to exhibit a true gap 
above Ep21. The width of this hard gap, also called the 
minigapAi, has been calculated as a purely quantum ef- 
fect by using random matrix theory (RMT) and its value 
scales with the Thouless energy, Et — h/2ra, where Td is 
the average (classical) dwell time a particle stays in the 
billiard between successive Andreev reflections 2 ^. 

Since the existence of this gap is expected to be related 
to the chaotic nature of the electronic motion, many at- 
tempts have been undertaken to explain this result in 
semiclassical terms22r— , however this appeared to be 
rather complicated. Indeed a traditional semiclassical 
treatment based on the so-called Bohr-Sommerfeld (BS) 
approximation yields only an exponential suppression of 
the DoS22r— . This apparent contradiction of this pre- 
diction with the RMT one was resolved quite early by 
Lodder and Nazarov 2 ^ who pointed out the existence 
of two different regimes. The characteristic time scale 
that governs the crossover between the two regimes is 
the Ehrenfest time te ~ |lnfi|, which is the time scale 
that separates the evolution of wave packets following es- 
sentially the classical dynamics from longer time scales 
dominated by wave interference. In particular it is the 
ratio t = te/ta, that has to be considered. 

In the universal regime, r = 0, chaos sets in sufficiently 
rapidly and RMT is valid leading to the appearance of 
the aforementioned Thouless gap 27 . Although the Thou- 
less energy Et is related to a purely classical quantity, 
namely the average dwell time, we stress that the appear- 
ance of the minigap is a quantum mechanical effect, and 
consequently the gap closes if a symmetry breaking mag- 
netic field is applied 3 ^. Similarly if two superconductors 
are attached to the Andreev billiard, the size of the gap 
will depend on the relative phase between the two super- 
conductors, with the gap vanishing for a 7r-j unction 3 ^. 

The deep classical limit is characterized by r — > oo, and 
in this regime the suppression of the DoS is exponential 
and well described by the BS approximation. The more 
interesting crossover regime of finite Ehrenfest time, and 
the conjectured Ehrenfest time gap dependence of Ref . l28l 
have been investigated by various mean a 12 ' 21 ' 36 " — . Due 
to the logarithmic nature of te, investigating numeri- 
cally the limit of large Ehrenfest time is rather difficult, 
but a clear signature of the gap's Ehrenfest time depen- 
dence has been obtained^— for r < 1. From an an- 
alytical point of view RMT is inapplicable in the finite 
te regime 1 ^, therefore new methods such as a stochastic 
method 38 using smooth disorder and sophisticated per- 
turbation methods that include diffraction effects 3 ^ have 
been used to tackle this problem. On the other hand a 
purely phenomcnological model, effective RMT, has been 
develope d 37 ' 44 and predicts a gap size scaling with the 
Ehrenfest energy Ee — H/2te- Recently Micklitz and 
Altland 4 ^, based on a refinement of the quasiclassical ap- 
proach and the Eilenberger equation, succeeded to show 



the existence of a gap of width ttEe oc 1/t in the limit 
of large t>1. 

Consequently a complete picture of all the available 
regimes was still missing until recently when we treated 
the DoS semiclassically 45 following the scattering ap- 
proach 4 ^. Starting for r = and going beyond the diago- 
nal approximation we used an energy-dependent general- 
ization of the work^ 7 - on the moments of the transmission 
eigenvalues. The calculation is based on the evaluation 
of correlation functions also appearing in the moments 
of the Wigner delay times^ 8 -. More importantly, the ef- 
fect of finite Ehrenfest time could be incorporated in this 
framework 49 leading to a microscopic confirmation of the 
te dependence of the gap predicted by effective RMT. 
Interestingly the transition between t = and r — oo 
is not smooth and a second gap at ttE-e was observed 
for intermediate t, providing us with certainly the most 
clear-cut signature of Ehrenfest time effects. 

In this paper we extend and detail the results obtained 
ini£. First we discuss Andreev billiards and their treat- 
ment using RMT and semiclassical techniques. For the 
DoS in the universal regime (t = 0) we first delve into 
the work of Refs. I47ll48l before using it to obtain the gen- 
erating function of the correlation functions which are 
employed to derive the DoS. This is done both in the 
absence and in the presence of a time reversal symme- 
try breaking magnetic field, and we also look at the case 
when the bulk superconducting gap and the excitation 
energy of the particle are comparable. 

We then treat Andreev billiards connected to two su- 
perconducting contacts with a phase difference <p. The 
gap is shown to shrink with increasing phase difference 
due to the the accumulation of a phase along the trajec- 
tories that connect the two superconductors. Finally the 
Ehrenfest regime will be discussed, especially the appear- 
ance of a second intermediate gap for a certain range of 
t. We will also show that this intermediate gap is very 
sensitive to the phase difference between the supercon- 
ductors. 



II. ANDREEV BILLIARDS 

Since the treatment of Andreev billiards was recently 
reviewed in Ref. HH we just recall some useful details here. 
In particular the chaotic Andreev billiard that we con- 
sider is treated within the scattering approach™ where 
the NS interface is modelled with the help of a ficti- 
tious ideal lead. This lead permits the contact between 
the normal metal cavity (with chaotic classical dynam- 
ics) and the semi-infinite superconductor as depicted in 

Fig. 

Using the continuity of the superconducting and nor- 
mal wave function, we can construct the scattering ma- 
trix of the whole system. Denoting the excitation energy 
of the electron above the Fermi energy E-p by E and as- 
suming that the lead supports N channels (transverse 
modes at the Fermi energy), the scattering matrix of the 



(a) (b) 

FIG. 1: (a) The Andreev billiard consists of a chaotic normal 
metal (N) cavity attached to a superconductor (S) via a lead, 
(b) At the NS interface between the normal metal and the 
superconductor electrons are retroreflected as holes. 



FIG. 2: (a) An Andreev billiard connected to two supercon- 
ductors (Si, S2) at phases ±</>/2 via leads carrying Ni and N2 
channels, all threaded by a perpendicular magnetic field B. 
(b) The semiclassical treatment involves classical trajectories 
retroreflected at the superconductors an arbitrary number of 
times. 



whole normal region can be written in a joint electron- 
hole basis and reads 



( S(E) 
^ S*{-E) 



(1) 



where S(E) is the unitary N x N scattering matrix of 
the electrons (and its complex conjugate S*(—E) that 
of the holes). As the electrons and holes remain un- 
coupled in the normal region the off-diagonal blocks are 
zero. Instead, electrons and holes couple at the NS in- 
terface through Andreev reflection^ where electrons are 
retroreflected as holes and vice versa, as in Fig. For 
energies E smaller than the bulk superconductor gap A 
there is no propagation into the superconductor and if we 
additionally assume A <C E-p we can encode the Andreev 
reflection in the matrix 



S A (E) = a(E) 



1 

1 



(2) 



where the actual setup treated is depicted in Fig. [3h.. 
It consists of a normal metal (N) connected to two su- 
perconductors (Si, S2) by narrow leads carrying N± and 
N2 channels. The superconductors' order parameters are 
considered to have phases ±0/2, with a total phase differ- 
ence 4>- Moreover a perpendicular magnetic field B was 
applied to the normal part. We note that although this 
figure (and Fig. [T^,) have spatial symmetry the treatment 
is actually for the case without such symmetry. 

As above, the limit A <C Ep was taken so that normal 
reflection at the NS interface can be neglected and the 
symmetric case in which both leads contain the same 
number, N/2, of channels was considere d 27 ' 35 . Finally it 
was also assumed that a — i, valid in the limit E, Et "C 
A -C E-p. For such a setup, the determinantal equation 
(|4]) becomes 



det l + S{E)e ict 'S*{-E)e^'t' =0 



(5) 



a(E) 



—1 arccos 



(«) 
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E 2 
A 2 "' 



(3) 



The retroreflection (of electrons as holes with the 
same channel index) is accompanied by the phase shift 
arccos (E/ A). In the limit of perfect Andreev reflection 
(E = 0) this phase shift reduces to tt/2. 

Below A the Andreev billiard has a discrete excitation 
spectrum at energies where det [1 — Sa(E)S^(E)] = 0, 
which can be simplified — to 



det [1 - a 2 (E)S(E)S*(-E)] = 0. 



(4) 



Finding the roots of this equation yields the typical den- 
sity of states of chaotic Andreev billiards. In the next two 
Sections we review the two main analytical frameworks 
that can be used to tackle this problem. 



A. Random matrix theory 

One powerful treatment uses random matrix theory. 
Such an approach was initially considered in Refs. I27ll35l 



where cf> is a diagonal matrix whose first N/2 elements 
are <fi/2 and the remaining N/2 elements — </>/2. We note 
that though we stick to the case of perfect coupling here, 
the effect of tunnel barriers was also included in Ref. [2?]. 

The first step is to rewrite the scattering problem in 
terms of a low energy effective Hamiltonian H 



H = 



H 



irXX 1 



-ttXX t -H* 



(6) 



where H is the MxM Hamiltonian of the isolated billiard 
and X an M x N coupling matrix. Eventually the limit 
M — > 00 is taken and to mimic a chaotic system the 
matrix H is replaced by a random matrix following the 
Pandey-Mehta distributior 
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P(H) cx exp 
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The parameter a measures the strength of the time- 
reversal symmetry breaking so we can investigate the 
crossover from the ensemble with time-reversal symme- 
try, the Gaussian orthogonal ensemble (GOE), to that 
without, the Gaussian unitary ensemble (GUE). It is re- 
lated to the magnetic flux $ through the two-dimensional 
billiard of area A and with Fermi velocity v-p by 



Ma 2 



hv-f, 



N 



27t£ , tv / X 



(8) 



Here c is a numerical constant of order unity depending 
only on the shape of the billiard. The critical flux is then 
defined via 



Ma 2 



N ( $ 



h /2ttE t . 

& $ c « - -r 1 Ai. 

e \ nvF 



(9) 

The density of states, divided for convenience by twice 
the mean density of states of the isolated billiard, can be 
written as 



d(e) = -ImW(e), 



(10) 



where W(e) is the trace of a block of the Green function 
of the effective Hamiltonian of the scattering system and 
for simplicity here we express the energy in units of the 
Thouless energy e = E/Et- This is averaged by inte- 
grating over ([7]) using diagrammatic methods 50 , which 
to leading order in inverse channel number 1/N leads to 
the expression^ 



W(e) 




1 + W 2 (e) + 



Vl + W 2 (e) 



(11) 

where /3 = cos (0/2) and b — (<f>/$ c ) with the critical 
magnetic flux <£> c for which the gap in the density of states 
closes (at = 0). Equation (TTTJ may also be rewritten 
as a sixth order polynomial and when substituting into 
([1TJ]) , we should take the solution that tends to 1 for large 
energies. In particular, when there is no phase difference 
between the two leads (4> = 0, or equivalently when we 
consider a single lead carrying N channels) and no mag- 
netic field in the cavity (Q/Q c = 0) the density of states 
is given by a solution of the cubic equation 

e 2 W 3 (e) + AeW 2 {e) + (4 + e 2 )W(e) + Ae = 0. (12) 



B. Semiclassical approach 

The second approach, and that which we pursue and 
detail in this paper, is to use the semiclassical approxima- 
tion to the scattering matrix which involves the classical 
trajectories that enter and leave the cavity^ 1 -. Using the 
general expression between the density of states and the 
scattering matrij*^, the density of states of an Andreev 
billiard reads^ikS 3 - 



1 d 

d(E) = d Im— In det [1 

7T 8E 



S A (E)S N (E)} , (13) 



where d — N/2ttEt is twice the mean density of states 
of the isolated billiard (around the Fermi energy) . Equa- 
tion (fT"3|) should be understood as an averaged quantity 
over a small range of the Fermi energy or slight variations 
of the billiard and for convergence reasons a small imag- 
inary part is included in the energy E. In the limit of 
perfect Andreev reflection a{E) w — i, see ([3]), and (fi"3|) 
reduces to 



d(E) = d 



— Im--— Tr — 

7T dE ^ m 

m — 1 



iS*(-E) 
iS(E) 



(14) 

Obviously only even terms in the sum have a non-zero 
trace, and setting n = 2m, dividing through by d and 
expressing the energy in units of the Thouless energy 
e = E/Et, this simplifies to2£ 



d(e) 



2Im H 



■1)" dC{e,n) 
n de 



(15) 



Equation (|15p involves the correlation functions of n scat- 
tering matrices 



C(e,n) = ^Tr 



S* 



—\s(— 

2t a ) \2r d 



(16) 



where we recall that the energy is measured relative to 
the Fermi energy and that E^ = h/2rd involves the aver- 
age dwell time Td. For chaotic systems^ the dwell time 
can be expressed as = T^/N in terms of the Heisen- 
berg time Th conjugate to the mean level spacing {2/d). 

At this point it is important to observe that nonzero 
values of e are necessary for the convergence of the ex- 
pansion of the logarithm in ([T3|) that led to ([T5|) . On the 
other hand, we are particularly interested in small values 
of e which puts (Ti"5l) on the edge of the radius of con- 
vergence, where it is highly oscillatory. The oscillatory 
behavior and a slow decay in n is a direct consequence 
of the unitarity of the scattering matrix at e — (in fact 
later it can also be shown that 8 | e =o = m )- Thus a 
truncation of (|15|) will differ markedly from the predicted 
RMT gap, which was the root of the difficulty of captur- 
ing the gap by previous semiclassical treatment s 3 ^ 33 ' 34 . 
In the present work we succeed in evaluating the entire 
sum and hence obtain results which are uniformly valid 
for all values of e. 

Calculating the density of states is then reduced to 
the seemingly more complicated task of evaluating cor- 
relation functions semiclassically for all n. Luckily the 
treatment of such functions has advanced rapidly in the 
last few year o 47 ' 48 ' 55 " — and we build on that solid basis. 
We also note that determining C(e,n) is a more gen- 
eral task than calculating the density of states. Since 
the Andreev reflection has already been encoded in the 
formalism before (TT51) . the treatment of the C(e,n) no 
longer depends on the presence or absence of the super- 
conducting material, but solely on the properties of the 
chaotic dynamics inside the normal metal billiard. 
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FIG. 3: (a) The original trajectory structure of the correlation function C(e, 2) where the incoming channels are drawn on the 
left, outgoing channels on the right, electrons as solid (blue) and holes as dashed (green) lines, (b) By collapsing the electron 
trajectories directly onto the hole trajectories we create a structure where the trajectories only differ in a small region called an 
encounter. Placed inside the Andreev billiard this diagram corresponds to Fig.[2p. The encounter can be slid into the incoming 
channels on the left (c) or the outgoing channels on the right (d) to create diagonal-type pairs. 



In the semiclassical approximation, the elements of the 
scattering matrix are given by^i 

S m (E) £ A z e^ E »\ (17) 

where the sum runs over all classical trajectories ( start- 
ing in channel i and ending in channel o. S(;(E) is the 
classical action of the trajectory £ at energy E above the 
Fermi energy and the amplitude Aq contains the stability 
of the trajectory as well as the Maslov phases 5 -^. After we 
substitute (fTT)) into (JT5J) and expand the action around 
the Fermi energy up to first order in e using 8Sq / dE = Tf 
where is the duration of the trajectory £, the correla- 
tion functions are given semiclassically by a sum over 2n 
trajectories 

xe 1 j . (18) 

The final trace in (|16j) means that we identify i n +i = ii 
and as the electron trajectories Q start at channel ij and 
end in channel Oj while the primed hole trajectories Q go 
backwards starting in channel Oj and ending in channel 
ij + i the trajectories fulfill a complete cycle, as in Figs. [3^ 
and|Ui,d,g. The channels i%, . . . , i n will be referred to as 
incoming channels, while 0\, . . . , o n will be called outgo- 
ing channels. This refers to the direction of the electron 
trajectories at the channels and not necessarily to which 
lead the channel finds itself in (when we have two leads 
as in Fig. (2). 

The actions in (|18p are taken at the Fermi energy and 
the resulting phase is given by the difference of the sum 
of the actions of the unprimed trajectories and the sum of 
the actions of the primed ones. In the semiclassical limit 
of h -> (c.f. the RMT limit of M oo) this phase 
oscillates widely leading to cancellations when the aver- 
aging is applied, unless this total action difference is of 
the order of h. The semiclassical treatment then involves 
finding sets of classical trajectories that can have such a 
small action difference and hence contribute consistently 
in the limit h — > 0. 



III. SEMICLASSICAL DIAGRAMS 

As an example we show the original trajectory struc- 
ture for n — 2 in Fig. [3Jt, where for convenience we draw 
the incoming channels on the left and the outgoing chan- 
nels on the right so that electrons travel to the right and 
holes to the left (c.f. the shot noise in Refs. [59j - [61T) . Of 
course the channels are really in the lead (Fig. QJi) or ei- 
ther lead (Fig. [2J and the trajectory stretches involve 
many bounces at the normal boundary of the cavity. 
We draw such topological sketches as the semiclassical 
methods were first developed for transpor t 47 ' 55 ! 57 where 
typically we have (complex conjugate transpose) in- 
stead of S* (complex conjugate) in (|16l) . restricted to the 
transmission subblocks, so that all the trajectories would 
travel to the right in our sketches. Without the mag- 
netic field, the billiard has time reversal symmetry and 
S is symmetric, but this difference plays a role when we 
turn the magnetic field on later. An even more impor- 
tant difference is that in our problem any channel can be 
in any lead. 

To obtain a small action difference, and a possible con- 
tribution in the semiclassical limit, the trajectories must 
be almost identical. This can be achieved for example 
by collapsing the electron trajectories directly onto the 
hole trajectories as in Fig. Inside the open circle, 
the holes still 'cross' while the electrons 'avoid cross- 
ing', but by bringing the electron trajectories arbitrarily 
close together the set of trajectories can have an arbi- 
trary small action difference. More accurately, the exis- 
tence of partner trajectories follows from the hyperbol- 
icity of the phase space dynamics. Namely, given two 
electron trajectories that come close (have an encounter) 
in the phase space, one uses the local stable and unsta- 
ble manifolds^— to find the coordinates through which 
hole trajectories arrive along one electron trajectory and 
leave along the other, exactly as in Fig.[3}D (and Fig. [2b) . 
These are the partner trajectories we pick for and £2 
when we evaluate C(e, 2) from (|18|) in the semiclassical 
approximation. As the encounter involves two electron 
trajectories it is called a 2-encounter. An encounter can 
happen anywhere along the length of a trajectory. In 
particular, it can happen at the very beginning or the 
very end of a trajectory, in which case it is actually hap- 
pening next to the lead, see Figs. |3b,d. This situation 
is important as it will give an additional contribution to 
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FIG. 4: (a) The original trajectory structure of the correlation function C(e, 4) where the incoming channels are drawn on the 
left, outgoing channels on the right, electrons as solid (blue) and holes as dashed (green) lines. (d,g) Equivalent 2D projections 
of the starting structure as the order is determined by moving along the closed cycle of electron and hole trajectories, (b) 
By pinching together the electron trajectories (pairwise here) we can create a structure which only differs in three small 
regions (encounters) and which can have a small action difference, (e) Projection of (b) also created by collapsing the electron 
trajectories in (g) directly onto the hole trajectories. (c,f) Sliding two of the encounters from (b) together (or originally pinching 
3 electron trajectories together) creates these diagrams. (h,i) Resulting rooted plane tree diagrams of (e,f) or (b,c) defining the 
top left as the first incoming channel (i.e. the channel ordering as depicted in (e,f)). 



that of an encounter happening in the body of the bil- 
liard. We will refer to this situation as an 'encounter 
entering the lead'. We note that if an encounter enters 
the lead the corresponding channels must coincide and we 
have diagonal-type pairs (i.e. the trajectories are coupled 
exactly pairwise) though it is worth bearing in mind that 
there is still a partial encounter happening near the lead 
as shown by the Ehrenfest time treatmen t 60 ' 65 . 

To give a more representative example, consider the 
structure of trajectories for n — 4. For visualization pur- 
poses in Fig. [4ji the original trajectories are arranged 
around a cylinder in the form of a cat's cradle. The in- 
coming and outgoing channels are ordered around the 
circles at either end although they could physically be 
anywhere. Projecting the structure into 2D we can draw 
it in several equivalent ways, for example as in Fig. 0ti 
or0g, and we must take care not to overcount such equiv- 
alent representations. We note that the ordering of the 
channels is uniquely defined by the closed cycle that the 
trajectories form. To create a small action difference, 
we can imagine pinching together the electron (and hole) 
strings in Fig.0^,. One possibility is to pinch two together 



in three places (making three 2-encounters) as in Fig. \^p. 
A possible representation in 2D is shown in Fig.013, which 
can also be created by collapsing the electron trajectories 
directly onto the hole trajectories in Fig. 0g. Note that 
the collapse of the diagram in Fig. 0J1 leads to a different 
structure with three 2-encounters. However in general it 
is not true that the different projections of the arrange- 
ment in Fig. 2^, are in a one-to-one correspondence with 
all possible diagrams. 

From Figs. HJ^e we can create another possibility by 
sliding two of the 2-encounters together to make a 3- 
encounter (or alternatively we could start by pinching 
three trajectories together in Fig. [4^, as well as an ad- 
ditional pair) as in Fig. Ht,f. Finally we could combine 
both to a single 4-encounter. Along with the possibilities 
where all the encounters are inside the system, we can 
progressively slide encounters into the leads, as we did 
for the n = 2 case in Fig. [3J creating, among others, the 
diagrams in Fig. [5] 

Finally, we mention that so far we were listing only 
'minimal' diagrams. One can add more encounters to 
the above diagrams but we will see later that such ar- 



FIG. 5: Further possibilities arise from moving encounters into the lead(s). Starting from Fig. [4J; we can slide the 2-encounter 
into the outgoing channels on the right (called 'o-touching', see text) to arrive at (a,d) or the 3-encounter into the incoming 
channels on the left (called 'i-touching') to obtain (b,e). Moving both encounters leads to (c,f), but moving both to the same 
side means first combining the 3- and 2-encounter in Fig. 3}j into a 4-encounter and is treated as such. 



rangements contribute at a higher order in the inverse 
number of channels and are therefore subdominant. The 
complete expansion in this small parameter is available 
only for small values of n, see Refs. [56ll57ll59l . 



A. Tree recursions 

To summarize the previous paragraph, the key task 
now is to generate all possible minimal encounter ar- 
rangements (see, for example, Ref. |H for the complete 
list of those with n = 3). This is a question that was an- 
swered in Ref. |4?] where the moments of the transmission 
amplitudes were considered. The pivotal step was to re- 
draw the diagrams as rooted plane trees and to show that 
there is a one-to-one relation between them (for the dia- 
grams that contribute at leading order in inverse channel 
number) . To redraw a diagram as a tree we start with a 
particular incoming channel i\ as the root (hence rooted 
trees) and place the remaining channels in order around 
an anticlockwise loop (hence plane). Moving along the 
trajectory £i we draw each stretch as a link and each en- 
counter as a node (open circle) until we reach o\. Then 
we move along £J back to its hrst encounter and continue 
along any new encounters to %2 and so on. For example, 
the tree corresponding to Figs. HjD,e is drawn in Fig. [4}r 
and that corresponding to Figs. 2J;,f is in Fig. @J. Note 
that marking the root only serves to eliminate overcount- 
ing and the final results do not depend on the particular 
choice of the root. 

A particularly important property of the trees is their 
amenability to recursive counting. The recursions be- 
hind our treatment of Andreev billiards were derived in 
Ref. H3 and we recall the main details here. First we can 
describe the encounters in a particular tree by a vector 
v whose elements vi count the number of Z-encounters in 



the tree (or diagram) ; this is often written as 2 V2 3" 3 
An Z-encounter is a vertex in the tree of degree 21 (i.e. 
connected to 21 links). The vertices of the tree that corre- 
spond to encounters will be called 'nodes', to distinguish 
them from the vertices of degree 1 which correspond to 
the incoming and outgoing channels and which will be 
called 'leaves'. The total number of nodes is V = Yli>i v i 
and the number of leaves is 2n where n is the order of 
the correlation function C(e, n) to which the trees con- 
tribute. Defining L = Yli>i^ v i^ we can express n as 
n = (L — V + 1). Note that the total number of links is 
L+n which can be seen as I links trailing each Z-encountcr 
plus another n from the incoming channels. For example, 
the 2 1 3 1 tree in Fig. @J has L — 5, V — 2 and contributes 
to the n — 4 correlation function. We always draw the 
tree with the leaves ordered i\, 0\, . . . , i n , o n in anticlock- 
wise direction. This fixes the layout of the tree in the 
plane, thus the name 'rooted plane trees—. 

From the start tree, we can also move some encounters 
into the lead(s) and it is easy to read off when this is pos- 
sible. If an Z-encounter (node of degree 21) is adjacent to 
exactly I leaves with label i it may 'z-touch' the lead, i.e. 
the electron trajectories have an encounter upon enter- 
ing the system and the corresponding incoming channels 
coincide. Likewise if a 2Z-node is adjacent to I o-leaves it 
may 'o-touch' the lead. For example, in Fig. @J the top 
node has degree 6, is adjacent to 3 ?-leaves (including the 
root) and can i-touch the lead as in Figs.[5}),e. The lower 
encounter can o-touch as in Figs.[S^,,d. In addition, both 
encounters can touch the lead to create Figs. [5b, f. 

Semiclassically, we add the contributions of all the pos- 
sible trajectory structures (or trees) and the contribution 
of each is made up by multiplying the contributions of its 
constituent parts (links, encounters and leaves). First we 
count the orders of the number of channels N. As men- 
tioned in Ref. (see also Sec. IIVI below) the multiplica- 
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FIG. 6: The tree shown in (a) is cut at its top node (of degree 6) such that the trees (b)-(f) are created. Note that to complete 
the five new trees we need to add an additional four new links and leaves and that the trees (c) and (e) in the even positions 
have the incoming and outgoing channels reversed. 



tive contribution of each encounter or leaf is of order N 
and each link gives a contribution of order 1/N. Together 
with the overall factor of 1/N, see equation (jTHJ) , the to- 
tal power of 1/N is 7, the cyclicity of the diagram. Since 
our diagrams must be connected, the smallest cyclicity is 
7 = if the diagram is a tree. The trees can be gener- 
ated recursively, since by cutting a tree at the top node 
of degree 21 (after the root) we obtain 21 — 1 subtrees, as 
illustrated in Fig. [5] 

To track the trees and their nodes, the generating func- 
tion F(x, Zi, z a ) was introduced^ where the powers of 

• Xi enumerate the number of l-encounters, 

• Zij enumerate the number of Z-cncountcrs that i- 
touch the lead, 

• z ^i enumerate the number of Z-encounters that o- 
touch the lead. 

Later we will assign values to these variables which will 
produce the correct semiclassical contributions of the 
trees. Note that the contributions of the links and leaves 
will be absorbed into the contributions of the nodes hence 
we do not directly enumerate the links in the generating 
function F. Inside F we want to add all the possible 
trees and for each have a multiplicative contribution of 
its nodes. For example, the tree in Fig. 01 and its relatives 
in Fig. [5] would contribute 

X3X2 + Zi,3%2 + X 3 Z 0i2 + Z ii3 Z . 2 = (X 3 + Zi,3) (X2 + Z 0>2 ) . 

(19) 

A technical difficulty is that the top node may (if there 
are no further nodes) be able to both i-touch and o-touch, 
but clearly not at the same time. An auxiliary generat- 
ing function / = fix, Zi, z D ) is thus introduced with the 



restriction that the top node is not allowed to i-touch 
the lead. We denote by 'empty' a tree which contains no 
encounter nodes (like Fig.|BJi). An empty tree is assigned 
the value 1 {i.e. /(0) = 1) to not affect the multiplica- 
tive factors. To obtain a recursion for / we separate the 
tree into its top node of degree 21 and 21 — 1 subtrees 
as in Fig. [51 As can be seen from the Figure, I of the 
new trees (in the odd positions from left to right) start 
with an incoming channel, while the remaining I — 1 even 
numbered subtrees start with an outgoing channel, and 
correspond to a tree with the i's and o's are reversed. For 
these we use the generating function / where the roles of 
the z variables corresponding to leaves of one type are 
switched so / = f(x, z a , Zi). The tree then has the con- 
tribution of the top node times that of all the subtrees 
giving xiff 1 - 1 . 

The top node may also o-touch the lead, but for this 
to happen all the odd-numbered subtrees must be empty 
(i.e. they must contain no further nodes and end directly 
in an outgoing channel). When this happens we just 
get the contribution of z Dt i times that of the / — 1 even 
subtrees: z Qt if . In total we have 



1=2 

and similarly 

00 

f = i+sr[ xl fif-i +Zitl f-i 
1=2 



(20) 



(21) 



For F we then reallow the top node to i-touch the lead 
which means that the even subtrees must be empty and 



9 



a contribution of , giving 



(22) 



z=i 



if we let Zi i = 1 (and also z a \ = 1 for symmetry). Pick- 
ing an o-leaf as the root instead of an i-leaf should lead 
to the same trees and contributions so F should be sym- 
metric upon swapping Zi with z a and / with /. These 
recursions enumerate all possible trees (which represent 
all diagrams at leading order in inverse channel number) 
and we now turn to evaluating their contributions to the 
correlation functions C (e, n) . 



IV. DENSITY OF STATES WITH A SINGLE 
LEAD 

To calculate the contribution of each diagram, 
Refs. 15514571 used the ergodicity of the classical motion 
to estimate how often the electron trajectories are likely 
to approach each other and have encounters. Combined 
with the sum rul o 55 i 67 to deal with the stability ampli- 
tudes, Ref. [56| showed that the semiclassical contribution 
can be written as a product of integrals over the dura- 
tions of the links and the stable and unstable separations 
of the stretches in each encounter. One ingredient is the 
survival probability that the electron trajectories remain 
inside the system (these are followed by the holes whose 
conditional survival probability is then 1) which classi- 
cally decays exponentially with their length and the de- 
cay rate l/r<j = N/Tn. A small but important effect is 
that the small size of the encounters means the trajec- 
tories are close enough to remain inside the system or 
escape (hit the lead) together so only one traversal of 
each encounter needs to be counted in the total survival 
probability 



exp 



N 

ty 



L+n 



^ ^ ~t~ ^ ^ ta ; 



(23) 



a = l 



where the tj are the durations of the (n + L) link stretches 
and t a the durations of the V encounters so that the 
exposure time t x is shorter than the total trajectory time 
(which includes / copies of each /-encounter). 

As reviewed in Ref. the integrals over the links and 
the encounters (with their action differences) lead to sim- 
ple diagrammatic rules whereby 

• each link provides a factor of Th/ [N (1 — ie)] , 

• each /-encounter inside the cavity provides a factor 
oi-N(l-Ue)/T^, 

with the (1 — i/e) deriving from the difference between 
the exposure time and the total trajectory time. Recall- 
ing the prefactor in (|18|) and that L is the total number 
of links in the encounters, it is clear that all the Heisen- 
berg times cancel. The channel number factor N~ 2n from 



these rules and the prefactor (with n = L — V + l) cancels 
with the sum over the channels in (IT51) as each of the 2n 
channels can be chosen from the N possible channels (to 
leading order). 

With this simplification, each link gives (1 — ie) , 
each encounter — (1 — iZe) and each leaf a factor of 1. To 
absorb the link contributions into those of the encounters 
(nodes) we recall that the number of links is n-\-^2 
where a labels the V different encounters. Therefore the 
total contribution factorizes as 



f 



n 



(l-ii. 



(24) 



Moving an /-encounter into the lead, as in Fig. [5] means 
losing that encounter, / links and combining / channels 
so we just remove that encounter from the product above 
(or give it a factor 1 instead). 



A. Generating function 

Putting these diagrammatic rules into the recursions 
in Sec. IIIIAl then simply means setting 



fl-i/e) 



xi 



z it i = z 0i i = l-f l -\ (25) 



(1 - ^) 



where we additionally include powers of f to track the 
order of the trees and later generate the semiclassical 
correlation functions. The total power of f of any tree is 
Sz>iC — ^) Vl = L — V = n—1. To get the required pref- 
actor of (1 — ie) _n in ([M]) we can then make the change 
of variable 



/ = p(l-ie), f = 



1 - ie' 

so that the recursion relation (1201) becomes 



S(l-ie) = l- 



1=2 



(26) 



(l-UeJ+jy-V -1 . (27) 

1=2 



and similarly for g. Using geometric sums (the first two 
terms are the / = 1 terms of the sums) this is 
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1 



1 ~ rgg (1 _ r gg) 1 - rg 



(28) 



We note that the since / is obtained from / by swapping 
Zi and z and in our substitution (f25|) Zj = z , the func- 
tions / and / are equal. Taking the numerator of the 
equation above and substituting g = g leads to 



1 



9 



rg- 



1 - ie 1 - ie 



-[g-l-ie] 



(29) 



To obtain the desired generating function of the semi- 
classical correlation functions we set F — G (1 — ie) in 
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([22]) , along with the other substitutions in ([25]) and (l26l) , 



G(e,r) 



1 — rg 



G(e,r)=J2r n - 1 C(e,n), (30) 



so that by expanding g and hence G in powers of r we 
obtain all the correlation functions G(e,n). This can be 
simplified by rearranging ([30]) and substituting into ([29]) 
to get the cubic for G directly 

r(r-l) 2 G 3 +r(3r+ie-3)G 2 + (3r+ie-l)G+l = 0. (31) 



when expanded in powers of r, agrees for a range of values 
of n with the expansion of (|33[) derived from the corre- 
lation functions obtained from (I3TT) . To show that ([33)) 
agrees with (|33p to all orders in r we use a differentia- 
tion algorithm to find an equation for the intermediate 
generating function 



I(e,r) 
I(e,r) 



1 <9G(e,r) _ d[rH(e,r)} 
i <9e <9r 



x -> r 



dC(e,n) 



de 



(36) 



B. Density of states 

The density of states of a chaotic Andreev billiard with 
one superconducting lead (|15|) can be rewritten as 



(32) 



a ^ (-l)«- 1 C7(e,n) 



where without the 1/n the sum would just be G(e, —1) in 
view of ([30]) . To obtain the 1 /n we can formally integrate 
to obtain a new generating function H(e,r), 



H{e ' r) = iik I c; "' r](lr - 



r™- 1 SG(e,n) 
in <9e 



so the density of states is given simply by 
d(e) = 1 - 2RcH(e,-l). 



(33) 



(34) 



To evaluate the sum in (|32p we now need to integrate 
the solutions of (|31l) with respect to r and differentiate 
with respect to e. Since G is an algebraic generating 
function, i.e. the solution of an algebraic equation, the 
derivative of G with respect to e is also an algebraic gen- 
erating function 68 . However, this is not generally true for 
integration, which can be seen from a simple example of 
/ = 1/x, which is a root of an algebraic equation, unlike 
the integral of /. Solving equation (|3"Tj) explicitly and 
integrating the result is also technically challenging, due 
to the complicated structure of the solutions of the cubic 
equations. Even if it were possible, this approach would 
fail in the presence of magnetic field, when G is a solution 
of a quintic equation, see Sec. IIVD1 or in the presence of 
a phase difference between two superconductors. 

The approach we took is to conjecture that H(e, r) 
is given by an algebraic equation, perform a computer- 
aided search over equations with polynomial coefficients 
and then prove the answer by differentiating appropri- 
ately. We found that 



{erf (I - r)H 3 + ier[r(ie - 2) + 2(1 - ie)]H 2 
+ [r(l-2ie)-(l-ie) 2 ]iJ+l = 0, 



both starting from (I3TT) and from (|35p and verifying that 
the two answers agree. 

The differentiation algorithm starts with the algebraic 
equation for a formal power series 77 in the variable x 
which satisfies an equation of the form 

$(x,ri) := p (x) + Pl (x)r] + . . . + p m (x)r] m = 0, (37) 

where po(x), • • • ,Pm{x) are some polynomials, not all of 
them zero. The aim is to find an equation satisfied by 
£ = d^/dx, of the form 



q {x) + qi (x)£ + . . . + q m (x)C = 0, 



(38) 



where qo(x), . . . , q m (x) are polynomials. Differentiating 
(|37|) implicitly yields 



d$(x, 77) / d$(x, rj) 



dx 



drj 



Pijhx) 
Q(i],x)' 



(39) 



(35) 



where P and Q are again polynomial. After substitut- 
ing this expression into the algebraic equation for £ and 
bringing everything to the common denominator we get 

q (x)Q m (x, 77) + Ql (x)P(x, ri)Q m -\x, 77) 

+ ... + q m (x)P m (x,r 1 ) = 0. (40) 

However, this equation should only be satisfied mod- 
ulo the polynomial &(x,rj). Namely, we use poly- 
nomial division and substitute P J (x 1 r))Q m ~ :1 (x, 17) = 
T(a;,77)$(a;,77) + Rj{x,rj) into ijlTj]). Using ([37]) we arrive 
at 

q (x)Ro(x,r)) +q 1 (x)R 1 (x,r)) + . . . + q m {x)R m {x,rf) = 0. 

(41) 

The polynomials Rj are of degree of m — 1 in rj. Treat- 
ing (j41[) as an identity with respect to rj we thus obtain 
m linear equations on the coefficients qj. Solving those 
we obtain qj as rational functions of x and multiplying 
them by their common denominator gives the algebraic 
equation for £. 

Performing this algorithm on G from pip , with x = ie, 
and on rH from l|35p. with x = r, leads to the same 
equation, given as (|A.1|) in the appendix, for the inter- 
mediate function defined in ([36)) and therefore proves the 
validity of the equation ([351) . Setting e = in ([35]) then 
= in as mentioned in Sec. Ill Bl 



shows that 



dC{e,n) 1 
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FIG. 7: (a) The density of states of a chaotic quantum dot coupled to a single superconductor at E <C A. (b) The density of 
states with a finite bulk superconducting gap A = 2Et (dashed line) and A = 8Et (solid line) compared to the previous case 
in (a) with A — > oo (dotted line). 



To compare the final result ([34]) with the RMT predic- 
tion we can substitute H(e, —1) = [— iW(e) + 1] /2 into 
([35]) . The density of states is then given in terms of W as 
d(e) = -ImH^(e). The equation for W simplifies to the 
RMT result (IT21) , and the density of states then reads^ 



d(e) 



e < 2 



s/3 



«+(«)- «-(«)] e >2(^ti 



5/2 
5/2 



(42) 

where Q±(e) = (8 - 36e 2 ± 3e%/3e 4 + 132e 2 - 48) V3 . 
This result is plotted in Fig. [7^, and shows the hard gap 
extending up to around 0.6-Et- 



C. Small bulk superconducting gap 

The calculation of the density of states above used the 
approximation that the energy was well below the bulk 
superconductor gap, E -C A or e -C 6 (for 6 = A/Et), 
so that the phase shift at each Andreev reflection was 
arccos(e/<5) ~ n/2. For higher energies or smaller super- 
conducting gaps, however, the density of states should be 
modified-^ to 



d{e) = 1 + Re 



2Im]T 



a(e) 2n C(e,n) 



(43) 

where a(e) = S/(e + iV S 2 — e 2 ) as in ([3]). When taking 
the energy derivative in the sum in (1431) we can split the 
result into two sums and hence two contributions to the 
density of states 



d(e) = l + 2Im^ 

n=l 

+ Rc 1 



a(e) 2n dC(e,n) 



VP 



de 

n=l 



(44) 



a(e) 2n C{e,n) 



Here the first term, which comes from applying the en- 
ergy derivative to C(e, n), gives an analogous contribu- 
tion to the case E -C A but with r — a 2 instead of — 1 
and involving H(e,a 2 ) from ([53")) and (pIS")) . The second 
term in (]44[) comes from the energy derivative of a 2n and 
can be written using G(e, a 2 ) from (|30f and (pTj) : 



d(e) = Re[l + 2a 2 J ff(e,a 2 )] 

-,=L=[l-|-2a 2 G( e ,a 2 )]. (45) 



Re 



The effect of a finite bulk superconducting gap on the 
hard gap in the density of states of the Andreev billiard 
is fairly small, for example as shown in Fig. 03 even for 
5 = A/Et = 2 the width just shrinks to around 0.5£t- 
For 5 — 2 the shape of the density of states is changed 
somewhat (less so for 6 = 8) and we can see just before 
e = 2 it vanishes again giving a second thin gap. This 
gap, and even the way we can separate the density of 
states into the two terms in (|45|) , foreshadows the effects 
of the Ehrenfest time (in Sec. IVI[) . For energies above 
the bulk superconducting gap (e > S) we see a thin sin- 
gular peak from the V 5 2 — e 2 which quickly tends to the 
density of states of an Andreev billiard with an infinite 
superconducting gap as the energy becomes larger. 



D. Magnetic field 

If a magnetic field is present, the time reversal symme- 
try is broken and we wish to treat this transition semi- 
classically as in Refs. Ifj4ll70l . Note that since for the lead- 
ing order diagrams each stretch is traversed in opposite 
directions by an electron and a hole we are effectively 
considering the same situation as for parametric correla- 
tions^^ 2 -. Either way, the idea behind the treatment is 
that the classically small magnetic field affects the classi- 
cal trajectories very little, but adds many essentially ran- 
dom small perturbations to the action. The sum of these 
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FIG. 8: The effect of a time reversal symmetry breaking mag- 
netic field on the density of states of a chaotic Andreev bil- 
liard with a single superconducting lead for 6 = (dotted 
line), 6 = 1/4 (solid line), 6=1 (dashed line) and 6 = 9/4 
(dashed dotted line). 



fluctuations is approximated using the central limit the- 
orem, and leads to an exponential damping so the links 
now provide a factor of Tn/N(l — ie + 6). The parameter 
6 is related to the magnetic field via b = (<!>/<I> c ) as in 
Sec. Ill Al For an ^-encounter however, as the stretches 
are correlated and affected by the magnetic field in the 
same way, the variance of the random fluctuations of all 
the stretches is I 2 that of a single stretch. Hence each en- 
counter now contributes N (l — Ue + l 2 b) /T^ and again 
the correlation inside the encounters leads to a small but 
important effect. 

Similarly to the treatment without the magnetic field 
above, we can put these contributions into the recursions 
m Sec. IfflAl by setting 



- (1-ik 



xi 



and 



l 2 b) 



(l-ie + 6) 1 



f l ~\ 



f = g(l-ie + b), 



1 



l.f l -\ 



(46) 



(47) 



The intermediate generating function is then given by the 
implicit equation 



„2„5 



-r'g" + (1 + ie + 6)r ,g 4 + (2 - ie - b)rg 6 
~ (2 + ie - b)rg 2 - (1 - ie + b)g + 1 = 0, (48) 

and the generating function G(e, 6, r) of the magnetic 
field dependent correlation functions G(e, 6, n), which is 
still connected to g via G — g/(l — rg), is given by 

r 2 (r- 1) 3 G 5 

+ (ier - ie + 5r 2 - lOr + 5 - br - b) r 2 G 4 
+ (3ier - ie + 10r 2 - 12r + 2 - 36r - 6) rG 3 
+ (3ie + lOr - 6 - 36)rG 2 

- (1 - 5r-ie + 6)G + 1 = 0. (49) 



Removing the magnetic field by setting b = reduces 
both of these equations (after factorizing) to the pre- 
vious results (|2"9"|) and (pTTj) . Next we again search 
for and verify an algebraic equation for H(e, 6, r) — 
l/(ir) J[dG(e, b, r)/de]dr, though the higher order makes 
this slightly more complicated, finding 

46V (r - 1) H 5 + 4br 3 [ie - 36 + r (2b - ie)] H 4 
+ r 2 [e 2 (1 - r) + 2ie6 (5 - 3r) - b (136 + 4) 
+ br (56 + 4) ] H s 

2 (ie - 36) (1 - ie + 6) (50) 



+ r((l-ie + 6) 2 + 46-l) 



H- 



(1 - ie + 6) 2 - r (1 - 2ie + 26) H + 1 = Q 



In order to check the agreement with the RMT result 
we substitute H(e,b,-1) = [-iW(e,b) + 1] /2 into (5DJ) . 
This leads to 

b 2 W 5 - 2beW A - (46 - 6 2 - e 2 ) W 3 + 2(2 - b)eW 2 
+ (4-46 + e 2 )W + 4e = 0, (51) 

which corresponds to the RMT result (fTTj) with no phase 
(0 = 0). The density of states calculated from this equa- 
tion is shown in Fig. [S] for different values of 6. The gap 
reduces for increasing 6, closes exactly at the critical flux 
(6 = 1) and the density of states becomes flat (at 1) as 
6 — > oo. 



V. DENSITY OF STATES WITH TWO LEADS 

Next we consider a classically chaotic quantum dot 
connected to two superconductors with Ni and N2 chan- 
nels respectively and a phase difference 0, as depicted 
in Fig. [2Ja,. The density of states, as in Sec. Ill Al and 
Refs. l35l)69L can then be reduced to equation (fl"5"j) but 
with 



C(e,(f),n) = -^Tr 



eh 
2r d 



(52) 

where <j) is again a diagonal matrix whose first N\ ele- 
ments from the first superconductor Si are (f>/2 and the 
remaining N2 elements from S2 are —<j>/2. Note that the 
case = corresponds to the previous case of a single 
superconductor with N = Ni + N2 channels. When we 
substitute the semiclassical approximation for the scat- 
tering matrix (|17l) into (|52|) , and especially if we write the 
scattering matrix in terms of its reflection and transmis- 
sion subblocks, the effect of the superconductors' phase 
difference becomes simple. Namely, each electron (un- 
primed) trajectory which starts in lead 1 and ends in 
lead 2 picks up the phase factor exp(— i(j>) while each un- 
primed trajectory going from lead 2 to lead 1 receives the 
factor exp(i0). Reflection trajectories which start and 



13 




FIG. 9: The paths may start and end in either of the two leads 
as shown. £4 as it travels from lead 1 to lead 2 obtains a phase 
factor exp(— i</>), £2 traveling back contributes exp(i0) while 
the others does not contribute any phase. The encounters 
are again marked by circles and Si and S2 denote the two 
superconducting leads at the corresponding superconducting 
phases ±0/2. This diagram is equivalent to the one in Fig. [If. 



end in the same lead have no additional phase factor, as 
depicted in Fig. [9] Since exchanging the leads gives the 
opposite phase, we expect the solution to be symmetric 
if we simultaneously exchange N\ with N 2 and change <f> 

to —<j). 

As these factors are multiplicative, we can equivalently 
say that each electron trajectory leaving superconductor 
1 or 2 picks up exp(— i0/2) or exp(i</>/2) while each one 
entering lead 1 or 2 picks up exp(i</>/2) or exp(— i(f>/2). 
To include these factors in our semiclassical diagrams, 
we can simply remember that in our tree recursions in 
Sec. lIII Al the channels we designated as 'incoming' chan- 
nels have electrons leaving them while electrons always 
enter the outgoing channels. Each incoming channel (in 
the original channel sum in (|18p) can still come from the 
N possible channels, but with the trajectory leaving it 
now provides the factor Ni exp(— i^/2) + N2 exp(i^/2). 
Similarly each outgoing channel now provides the com- 
plex conjugate of this factor. Recalling the power of 
N~ 2n coming from the links and encounters, we can up- 
date the contribution of each diagram or tree (|24|) to 



iVie~~ + N 2 e~ 



') ( 



iVie~ + N 2 e~~ 



N 2n _ ie y 
(1 - U a e) 



S -i')' 



(53) 



However, moving an Z-encounter into lead 1 means com- 
bining I incoming channels, I links and the encounter 
itself. These combined incoming channels, with I elec- 
tron trajectories leaving, will now only give the factor 
Ni exp(— il(f>/2) + N 2 exp(iZ</>/2) where the important dif- 
ference is that I is inside the exponents. We therefore 
make the replacement 



when we move the encounter into the outgoing leads we 
take the complex conjugate of ([53j). 

To mimic these effects in the semiclassical recursions 
we can set 



r 



(1-ie)' 
(Nie~% + N 2 e% 



ZoJ 



f = 9 



( 



N 



(55) 



TVie" 



N 2 e~ 



N(3 l 
(Nxe^ + N 2 e~ 



N(/3*y 



f l -\ 



r l -\ 



P/3* ' 



(I -iey 



(56) 
(57) 



in Sec. IIII Al Including these substitutions in the recur- 
sion relation (|20p and summing we obtain 



if/: 



'9 



1 



(3(3* 



rgg 



((3(3* - rggf 
N 2 1 



N 



(3*e~ 



rg 



N 



(58) 



rg 



and a similar equation from (|2 1[) . The generating func- 
tion of the correlation functions C(e,(f>,n) is then given 
from ([22]) by 



G 



Ni 



N (3e'-f 



+ 



No 



fJ 



rg 



N (3e^ 



(59) 



rg 



Returning to (f58|) and multiplying through by g, we 
can see that the first two terms are symmetric in g and 
g. Combining the other two and taking the difference 
from the corresponding equation for g we have 



9 


\/3*f - rg 




(j3*e 2 — rg 


(^(3*e 2 — rg^j 



rg 



y(3e 2 — rg^j (j3e 2 — rgj 



(60) 



The resulting quadratic equation, when substituted back 
into (l58l) leads to a sixth order equation for g. Note that 
the right hand side of (|6"0")) is (recalling (|55|) and that 
N% + N 2 = N) the same as (|5"9")l so it is clear that G 
satisfies the required symmetry upon swapping the leads 
(i.e. swapping Ni with N 2 and <p with —(/))■ 



Aic~— + N 2 e^~ 
N 



(54) 



as well as removing the encounter from (I53|) . Similarly 



A. Equal leads 

To make the equations more manageable we focus for 
now on the simpler case in which the leads have equal 
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FIG. 10: The density of states of a chaotic quantum dot cou- 
pled to two superconductors with the same numbers of chan- 
nels and phase differences (dotted line), 57r/6 (solid line), 
21tt/22 (dashed line) and 123tt/124 (dashed dotted line). 
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FIG. 11: Magnetic field dependence of the density of states 
of a chaotic Andreev billiard with phase difference <j> — 5ir/6 
for 6 = (dotted line), b = 0.1024 (solid line), b = 0.4096 
(dashed line) and b = 1 (dashed dotted line). 



size and TVi = N 2 = N/2. Then (3 = cos(0/2) is real and 
we can see from ([50]) or = z Q that g = g is a solution. 
Putting this simplification into (|58p we can obtain the 
following quartic 

r 2 g 4 -r(l + r + ier)g 3 + 2ie/3 2 rg 2 + (l-ie + r)/3 2 g-l3 4 = 0. 

(61) 

We may also find an algebraic equation of fourth order 
for G if we solve ([5TJ1) for g and substitute the solution 



f3 2r/3G + f3- ^ f3 2 + 4rG (I + rG) (/3 2 
2 



1) 



r(l+rG) 



(62) 



into (IBTj) . Note that we take the negative square root to 
agree with the previous result when the phase is (i.e. 
(3=1) though this sign does not affect the equation one 
finally finds for G. After the fourth order equation for 
G has been found we can again search for and verify an 
equation for H(e, (j), r) = l/(ir) J(dG(e, <j>, r)/de)dr, 

eV [1 - 2r (2/3 2 - l) + r 2 } H 4 
+ ier 2 [2 - 3ie - 4r (1 - ie) (2f3 2 - l) 
+ r 2 (2 - ie)] H 3 

- r [1 - 4ie - 3e 2 - 2r (l - 3ie - e 2 ) (2/3 2 - l) 

+ r 2 (1 - 2ie)] i? 2 

- (1 - ie) 2 - 2r (1 - ie) (2/3 2 - l) + r 2 j i? 

+ /3 2 = 0. (63) 

In order to see the agreement of our result with the 
RMT prediction we again substitute H(e,<j>,—1) = 
[-iW(e,<£) + l]/2 such that d(e) = -lmW(e,4>). If we 
do so we find 



e 2 p 2 W A + 4eP 2 W 3 + [Af3 2 
+ AeP 2 W - e 2 + e 2 /? 2 = 



2e 2 / 3 2 )W 2 



(64) 



which corresponds to (|11[) for zero magnetic field. More- 
over, if the phase difference is zero (and /3 = 1), we can 
take out the factor W and recover IT2"]) . 

Solving this equation yields the density of states. If we 
insert different values for the phase cj> one finds that the 
hard gap in the density of states decreases with increasing 
phase difference while the density of states has a peak at 
the end of the gap which increases and becomes sharper 
with increasing phase. Finally when the phase difference 
is equal to ir the gap closes and the peak vanishes so the 
density of states becomes identical to 1. This can all be 
seen in Fig. [TU1 



B. Magnetic field. 

In the presence of a magnetic field, we again have to 
change the diagrammatic rules as in Sec. IIVDI Doing 
the calculation above with these modified diagrammatic 
rules leads to a sixth order equation for g: 

r 3 g e - r 2 [i + r (i + i e + b)} g 5 - r 2 [3 2 (1 - 2ie - 2b) g 4 
+ r/3 2 [2 - ie - b + r (2 + ie - 6)] g 3 
- r(3 4 (1 + 2ie - 2b) g 2 

-/3 4 {l + r-ie + b)g + p 6 = 0. (65) 

The relation (1591) between G and g remains unchanged 
and therefore we may find a sixth order equation for G. 
We find the corresponding if, which is recorded as (|A.2[) 
in the appendix, using a computer search over sixth or- 
der equations with polynomial (in e, (j), b and r) coeffi- 
cients whose expansion in r Q33p matches the correlation 
functions calculated by expanding G. We note that for 
this order polynomial it was not feasible (in terms of 
computational time and memory) to solve the equations 
resulting from the differentiation algorithm described in 



15 




d(e) 



0.8- 
0.6 
0.4 
0.2- 








0.5 



1.5 



(c) 



FIG. 12: Phase dependence of the density of states of a chaotic Andreev billiard with phase difference — (dotted line), 
4> — 7r/2 (solid line), </> — 5ir/6 (dashed line) and <f> = 2l7r/22 (dashed dotted line), (a) At magnetic field b = 0.1024, (b) at 
b = 0.4096 and (c) at b = 1. 



Sec. IIVBI and to find the intermediate generating func- 
tion / in all generality. However, we succeeded in finding 
a polynomial equation for / that was satisfied by the 
derivatives of both rH and G for a large number of nu- 
merical values of the parameters (e, 0, b). For each pa- 
rameter involved, the number of the values checked was 
larger than the maximum degree of the parameter in the 
conjectured equation. While we cannot rule out the pos- 
sibility that the true equation for / has a higher order, 
given the large number of numerical values checked this 
is highly unlikely. 

From H we obtain the equation for W(e, 0, b), 

b 2 (3 2 W 6 - 2ebfi 2 W 5 + (2b 2 P 2 + e 2 /3 2 - Ab/3 2 - b 2 ) W 4 
+ 2 (eb + 2ef3 2 - 2eb(3 2 ) W 3 



+ (A/3 2 - b 2 - e 2 - Ab/3 2 + b 2 /3 2 + 2e 2 (i 2 ) W : 
+ 2(eb + 2ef3 2 - eb/3 2 ) W - e 2 + e 2 /3 2 = 0, 



(66) 

(ED 



which corresponds exactly to the full RMT result 
expanded. 

As an example, the magnetic field dependence of the 
density of states is shown at the phase difference of 57r/6 
in Fig. [TT] As the magnetic field is increased one finds a 
reduction of the gap and the peak appearing for a phase 
difference > vanishes again. Moreover the higher the 
phase difference, the lower the magnetic field needed to 
close the gap. While for = the gap closes at b = 1 
in the case of a phase difference of 57r/6 one needs b ps 
0.4096 and for = 2l7r/22 a magnetic field corresponding 
to b m 0.1024 closes the gap. In particular the critical 
magnetic field for which the gap closes is given by 35 



b B = 



2 cos (0/2) 
1 + cos (0/2)' 



(67) 



For ever increasing magnetic field the density of states 
approaches 1 and we can see that a higher phase dif- 
ference causes a faster convergence to this limit. Some 
examples are plotted in Fig. [T2"land there we see that for 
b = 1 the curve for = 2l7r/22 is nearly constant. 



C. Unequal leads 

Removing the restriction that the leads have equal size 
we return to a sixth order polynomial for g and G when 
substituting into (|55|) and then (|5^|) . Expanding 
G as a power series in r via G = r n ~ 1 C(e, 0, n) now 
gives three starting values for C(e, 0, 1) and we choose the 
one that coincides with the result from the semiclassical 
diagrams, namely /3/3*/ (1 — ie). Choosing the variable y 
to represent the relative difference in the lead sizes 



V 



N x -N 2 
N ' 



j3 = cos 



vy sin 



(68) 



leads to a particularly compact solution, and as be- 
fore, we can go through our roundabout route of finding 
the generating function of interest H(e,<p,y,r), which is 
recorded as (IA.3I) in the appendix. Although it also was 
not possible to verify (other than at a large number of 
parameter values) this sixth order equation, from it we 
can obtain the polynomial satisfied by W(e, 0, y): 

e 2 (3 2 W 4 + 4e$ 2 W 3 + (i$ 2 - e 2 + 2e 2 /3 2 ) W 2 



+ Aef3 2 W - e 2 + e 2 /3 2 
0, 



(2 + eJ-F) : 



4eV (l-/? 2 ) 



(69) 



where we have defined $ = cos(0/2) as the real part of 
j3 (which is equal to {3 when the leads have equal size) 
and the evenness in y follows from the symmetry under 
swapping the leads and to —0. The term in the square 
brackets is simply (|64[) and so we recover the result with 
equal leads when y — 0. Likewise we can check that when 
we only have a single lead (y = ±1) we recover a factor 
corresponding to (|T2l so that the phase, as expected, no 
longer plays a role. From this equation we can plot the 
density of states as in Fig. [T3] and see how the difference 
in lead sizes y interpolates between the result with equal 
leads above and the density of states with a single lead in 
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FIG. 13: Dependence of the density of states of an Andreev billiard on the difference y = (Ni — N2) /N in size of the leads 
with y — (dashed dotted line), y = 4/5 (dashed line), y = \/24/5 (solid line) and y = 1 (dotted line), (a) At phase difference 
cj> = 2tt/3, (b) at <f> = 5tt/6 and (c) at 4> = 21tt/22. 
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FIG. 14: Critical value of the difference in the lead sizes y 
as a function of the phase difference <j> between the two leads 
above which a second gap appears in the density of states. 



(|42j) . Note in particular that the peak in the density of 
states as the phase difference nears it vanishes slowly as y 
approaches ±1 so that we can see a second gap appear in 
the density of states for leads differing distinctly in chan- 
nel numbers (for example, see the solid line in Figs. [T3"b 
and c). Numerically we can extract the critical value of y 
for each phase difference (j> above which we see a second 
gap. We plot this in Fig. [14] where we see that the sec- 
ond gap only appears for particularly unequal leads and 
at reasonable phase differences. 



VI. EHRENFEST TIME DEPENDENCE 

So far we have been looking at the regime where 
the Ehrenfest time te ~ |hi/i|, the time below which 
wave packets propagate essentially classically (and above 
which wave interference dominates), is small compared to 
the dwell time t<i, the typical time the trajectories spend 
inside the scattering region. This is the same limit de- 
scribed by RMT and we have seen the agreement between 
semiclassics and RMT in Sees. II Vl and [Vl above. Moving 



away from this limit we can treat the typical effect of the 
Ehrenfest time on the correlation functions C(e, n), for 
now for the simplest case of a single lead and no mag- 
netic field. To contribute in the semiclassical limit, the 
correlated trajectories should have an action difference 
of the order of h which in turn means that the encoun- 
ters have a duration of the order of the Ehrenfest time. 
Increasing this relative to the dwell time, or increasing 
the ratio r = te/t^, then increases the possibility that 
all the trajectories travel together for their whole length 
in a correlated band. Likewise the probability of forming 
the diagrams (as in Fig.0} considered before reduces. All 
told, the Ehrenfest time dependence^ leads to the simple 
replacement 

C(e, r, n) = C(e, ^e^ 1 "™^ + . (70) 

1 — me 

This replacement leaves the n = 1 term unchan ged and 
had previously been shown for n = 2 in Ref. |6(J and 
n = 3 in Ref. 139 . The exponential growth of differences 
between trajectories due to the chaotic motion means 
that we just add the first term from the previous dia- 
grams with encounters in (|70|) to the second term from 
the bands as their opposing length restrictions lead to a 
negligible overlap. In fact this separation into two terms 
was show n 73 i 74 to be a direct consequence of the splitting 
of the classical phase space into two virtually independent 
subsystems. 

We leave the technical demonstration of (1701) to Ref. 
but the result follows by treating the diagrams consid- 
ered before, which are created by sliding encounters to- 
gether or into the lead (like the process depicted in Figs. [4] 
and [5| , as part of a continuous deformation of a single 
diagram. With a suitable partition of this family one can 
see that each set has the same te dependence and hence 
that ([70]) holds for all n. It is clear that in the limit 
r = ([70)) reduces to the previous (and hence RMT) 
results while in the opposite limit, r = 00, substituting 
(|70p into (|32|) and performing a Poisson summation we 
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FIG. 15: (a) Density of states for r = r^/rd = 2 (solid line), along with the BS (dashed) limit r — > oo and the RMT (dotted) 
limit r = 0, showing a second gap just below er = n. (b) Ehrenfest time related 27r/r-periodic oscillations in the density of 
states after subtracting the BS curve. 
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FIG. 16: (a) Width (and end point) of the first gap 
width of the second gap as a function of r. 



obtain the Bohr-Sommerfeld (BS)^ result 

2 cosh(7r/e) 



and (b) 



d BS (e) = (- 



sinh 2 (7r/e) 



(71) 



This result was previously found semiclassically by 
Ref. HO and corresponds to the classical limit of bands 
of correlated trajectories. 

For arbitrary Ehrenfest time dependence we simply 
substitute the two terms in dTOj) into ([32]) . With the 
second term we include 1 — (1 + r)e _r from the con- 
stant term (this turns out to simplify the expressions) 
and again perform a Poisson summation to obtain 



da(e,T) = l-(l + r)e- T 

" (-1)" d 
n de 



(72) 



+ 21m 

^Bs(e) 
— exp 



n=l 



-(1— ine)r 



1 — ine 



2irk\ 



d B s{e) + 



2fc(7r/e) 2 
sinh(7r/e) 



where k — [{er + 7r)/(27r)J involves the floor function, 
and we see that this function is zero for er < n. 



Of course the first term in (|70l) also contributes and 
when we substitute into (j32j) we obtain two further terms 
from the energy differential. These however may be writ- 
ten, using our semiclassical generating functions, as 

di(e,T) = e" T [1 - 2Re e i£T ff(e, -e ieT )] (73) 
+ tc- t [1 - 2Re e ier G(e, -e ier )] . 

Because G and H are given by cubic equations, we can 
write this result explicitly as 



di(e,r) 



v^e 



6e 



-Re[Q+(c,T)-Q_(e,T)] (74) 
-Re[P+(e,T)-P-(e,T)}, 



where 
Q±(e,r) 



8 _24e(l r cos(er))_ 24e2 
sin(eT) 

24e 2 (1 - cos(er)) 6e 3 (1 - cos(er)) 

sin 2 (eT) sin(er) 
2e 3 (2-3 cos(er) + cos 3 (er)) 



+ 



± 



sin 3 (eT) 
6eV3D (1 - cos(er)) 



sin (er) 



36e 



9e 2 sin(er) 



(l + cos(er)) 2 (l + cos(er)) 3 



(75) 



(76) 



+ 



± 



3V3D 



(l + cos(er)) 3 (l + cos(er)) 2 



These all involve the same discriminant D and so the 
differences in (T74l are only real (and hence d±(e, r) itself 
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is non-zero) when 

D(e,r) = e 4 -8e 3 sin(er)+4e 2 [5 + 6cos(er)] 

+ 24e sin(er) - 8 [1 + cos(er)] , (77) 

is positive. Recalling that the second contribution is zero 
up to er = 7r, the complete density of states is therefore 
zero up to the first root of D(e,r). The width of this 
gap is then solely determined by the contribution from 
quantum interference terms given by the trajectories with 
encounters. The hard gap up to the first root shrinks as r 
increases (see Fig. 116a ) and when taking the limit r — > oo 
while keeping the product er constant ([77]) reduces to 
—8 [1 + cos(er)] which has its first root at er = n. The 
gap then approaches E = ttE-^ for r 3> 1 where E^ — 
H/2te is the Ehrenfest energy. So one indeed observes a 
hard gap up to ttE-e in the limit r — > oo at fixed er in 
agreement with the quasiclassical result of Ref. HO- 

Alongside this reduction in size of the first gap, which 
was predicted by effective RMT 13 , when r > 0.916 the 
discriminant (f77|) has additional roots. Between the sec- 
ond and third root D(e, r) is also negative and a second 
gap appears. As r increases the roots spread apart so the 
gap widens. For example, the complete density of states 
for r = 2 is shown in Fig. ITSa along with the oscillatory 
behavior visible at larger energies (with period 2tt/t) in 
Fig. 115b. There the second gap is clearly visible and only 
ends when the second contribution c?2(e, r) becomes non- 
zero at er = it. In fact for r > ir/2 the third root of 
D(e,r) is beyond er = tt so the second gap is cut short 
by the jump in the contribution c?2(e,r). Since the sec- 
ond root also increases with increasing r the gap shrinks 
again, as can be seen in Fig. [T5b . 

To illustrate this behavior further, the density of states 
is shown for different values of r in Fig. [T7J One can see 
first the formation and then the shrinking of the second 
gap. As can be seen in the inset of Fig. ITTb the second 
gap persists even for large values of r and the size of the 
first hard gap converges slowly to er = 7r. The plot for 



r = 20 also shows how the density of states converges to 
the BS result. 



A. Effective RMT 

As mentioned above, the shrinking of the first gap has 
been predicted by effective RMT where the effect of the 
Ehrenfest time is mimicked 3 -^ by reducing the number 
of channels in the lead by a factor e T (to correspond to 
the part of classical phase space where the trajectories 
are longer than the Ehrenfest time) and to multiply the 
scattering matrix by the phase e leT ^ 2 to represent the 
energy dependence of the lead. The random matrix dia- 
grammatic averaging leads to the set of equation s 13 ' 44 

TF 2 + 1 = Wl (78) 
PF + W^sinw = — {W 2 + cosu + W sinu) , 

where u — er/2 and the density of states is given by (for 
u < it/2) 

d(e,r) = ~e- T lm(w —W 2 ). (79) 

V cos u / 

The equations in ([75} can be simplified to obtain a cu- 
bic for W (and W2) and in this form we can compare 
with our semiclassical results. In fact, making the sub- 
stitution H = [iW — l]/2r and setting r = — exp(ier) 
to get the first part in (fT5|) in the form of the first 
term in (|79[) we obtain exactly the same polynomial 
and hence agreement. Likewise when we substitute 
G = — [1UW2/ cos u + t] /2rr we obtain the same polyno- 
mial for the second part, albeit with the real offset tan it 
which does not affect the density of states. 

Of course this agreement provides semiclassical sup- 
port for the phenomenological approach of effective RMT 
as well as showing that (IT9"|) is valid for u beyond ir/2. In 
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FIG. 18: Density of states for r = 2 (solid line) along with the r = (dotted) and r = oo (dashed) limits for a chaotic Andreev 
billiard with phase difference (a) <j> — tv/18, (b) <j> = 57r/6 and (c) = 7r. 



principle then the second gap could also have been found 
using effective RMT, but of course effective RMT misses 
the second contribution to the density of states ^(e, r). 
This contribution turns out to be straightforward to ob- 
tain semiclassically, and can be compared to the bands 
treated in Ref. l40l 



B. Two superconducting leads 

If we include the effect of a symmetry breaking mag- 
netic field then, because of the way this affects the con- 
tribution of different sized encounters (as described in 
Sec. lIVD]) . such a simple replacement as in (|70p no longer 
holds. This situation has however been treated using ef- 
fective RMT— allowing them to also determine how the 
critical magnetic field (at which the gap closes) depends 
on the Ehrenfest time. Without the simple replacement 
of the magnetic field dependent correlations functions we 
are currently unable to confirm this result semiclassically. 
But if we look at two superconducting leads (for simplic- 
ity of equal size) at different phase this effect can be 
included in the channel sum and treated as above (the 
effective RMT result can be found by a simple modifica- 
tion of the treatment in Ref. l44l) . Important to remember 
is that the second part (of ([70)) ) corresponds to bands of 
trajectories that are correlated for their whole length and 
so they all start and end together (in the same leads). 
Therefore the second contribution has to be multiplied 
by [1 + cos(n0)]/2 leading to 

C(e,cj),T,n) = C(e,0,n)e- (1 - ine)r (80) 
1 + cos(n0) 1 - e -V~™)T 



+ 



1 — ine 



The first part of the density of states for non zero Ehren- 
fest time then remains as in f|T3[) but with G(e, r) and 
H(e, r) replaced by G(e, 0, r) and H(e, 0, r), respectively. 
The second contribution in this case however may be 
written as the average of the = contribution and 



a contribution with the full phase difference 0, 

daM,T) = i[4(e,O,r)+4(e,0,r)]. (81) 

Here (^(e, 0, r) may be again written as the sum of the 
t = oo result 

7T 



4 (1) M,r) 



(7T + 27rfcl 



2e 2 sinn 2 (tt/c) 



(82) 



) cosh 



7T — 27Tfcl 



+ (tt — 2nki + 0) cosh 



7r + 2irki 



and some correction 



y(2) 



2e 2 simi (tt/c) 



(83) 



) sinh ^ — 



| Trcosh (-) + (2ttA; 2 
7T cosh ^— J + (27rfc3 + cf>) sinh 



with ki = L(tt + cj>) /(2tt)J , k 2 = [{er + vr + 4>) /(2tt)J and 
^3 = L( eT + n — 4>) /(2tt)J • Since the k% and 4> only occur 
in the combinations 27rfci — <fi, 2nk2 — 4> and 27r/c3 + it is 
obvious that these contributions have oscillations in the 
phase <fi with period 27r. It can also be easily seen that 
for cf> — the previous result for the density of states in 
the Ehrenfest regime is reproduced. 

With |0| < 7r we have fci = k 2 = k% = for et < 
tt — |0| . Therefore one again sees that d 2 — as long 

as er < 7r — \<f>\. The first part equals the Bohr- 
Sommerfeld result (|7ip , so in the limit r = oo this result 
is reproduced again. The oscillations in e seen in the 
= case which have a period of 27r/r can still be seen 
due to the fact that the = result enters d 2 (e,<j),T) 
even if ^ 0. However one gets additional (but smaller) 
steps at energies satisfying e = [(2m — l)7r =F 0]/ t f° r 
integer m. 

We plot the density of states for r = 2, along with the 
t = and t = oo limits in Fig. [IS] for different values of 
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FIG. 19: Density of states for r = 1/2 (dotted line), r = 1 (dashed) and r = 2 (solid) showing the phase dependent jumps for 
phase difference (a) <f> — tt/18 and (b) <f> = 57r/6. 



the phase difference. We can see that as the phase differ- 
ence increases the second intermediate gap (c.f. Fig. ITSk) 
shrinks quickly. The reason for this shrinking is twofold: 
On the one hand, the gap in the RMT-like contribution 
shrinks, and on the other hand, the second contribution 
is zero only up to er = n — \<p\. Moreover if — > 7r the 
modified correlation function tends to zero so the density 
of states converges to (1 + r)e~ r + ^(e, r). For a finer 
look at the Ehrenfest time dependence and the phase de- 
pendent jumps we plot the density of states for r = 1/2, 
1 and 2 for phases = 7r/18 and 5ir/6 in Fig. [121 



VII. CONCLUSIONS 

From the semiclassical treatment of the density of 
states of chaotic Andreev billiards we have seen how fine 
correlations between ever larger sets of classical trajec- 
tories lead to the interference effects which cause a hard 
gap in the density of states. This treatment (c.f. the 
reservations in Ref. 38) builds on the recent advances in 
identifying 55 , codifying 5 ^ 5 - 7 - and generating^ 7 - the semi- 
classical contributions, and, because of the slow conver- 
gence of the expansion for the density of states in (fl~5|) . 
relies on the ability to treat correlations between n tra- 
jectories for essentially all n. The correlations between 
these trajectories, encoded in encounter regions where 
they differ slightly, are represented by simple (tree) dia- 
grams. These diagrams are related to those that appear 
for the conductance 5 ^ say where for increasing n they 
cause ever decreasing (in inverse channel number) cor- 
rections; here though they all contribute with roughly 
the same (slowly decreasing) importance. It is because 
we need to treat all orders that Andreev billiards are so 
interesting and the resultant effects so large. 

Along with obtaining the minigap, found by RMT 27 , 
for a billiard with a single lead, we could also obtain the 
full result for the density of states of an Andreev billiard 



with two superconducting leads at phase difference </>, 
treated using RMT in Ref. [H. The semiclassical paths 
that connect the two leads accumulate phases e ±I( ^ and 
cause the gap to shrink with increasing phase difference. 
It was also possible to treat the effect of a time rever- 
sal symmetry breaking magnetic field fo, considered with 
RMT in Ref. |35|, which makes the formation of the clas- 
sical trajectory sets, traversed in opposite directions by 
an electron and a hole, less likely. This in turn leads to 
a reduction of the minigap and a smoothing of the den- 
sity of states, especially for large phase differences 4>. We 
have found that in the limits <j) — > tt and b — > oo quantum 
effects vanish and the density of states becomes identical 
to the density of states of the isolated billiard. 

The agreement shown here between the semiclassical 
and the RMT results may lead one to wonder about the 
deeper connections between the two methods. Indeed the 
diagrammatic methods^ used in RMT, which also use re- 
cursion relations over planar diagrams, bear some resem- 
blance to the tree recursions here, but unfortunately any 
correspondence between the two is somewhat obscured 
by the transformation from the generating function G to 
the generating function H. It is also worth pointing out 
that our semiclassical treatment (with its inherent semi- 
classical limits) gives us access to the typical and univer- 
sal density of states of chaotic systems. However there 
can be further effects that change the shape of the den- 
sity of states, for example scarring 7 - 5 -, classical Lyapunov 
exponent fluctuations^ and disorder—. 

Of course all our results in this paper (and the RMT 
one o 27 ' 35 ) are only valid to leading order in inverse chan- 
nel number. With the formalism shown in this paper, to 
go to subleading order we only require a way of generat- 
ing the possible semiclassical diagrams. The contribution 
of eac h 56 ' 57 and how they affect the density of states is 
known in principle, but the key problem is that the struc- 
ture we used here breaks down, namely that in the tree 
recursions when we cut a rooted plane tree at a node we 
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created additional rooted plane freest. How to treat the 
possible diagrams which include closed loops etc, though 
generated for n = 1 in Ref. and n — 2 in Ref. H3 
by cutting open closed periodic orbits, remains unclear. 
However the treatment for n = 1 and n = 2 makes it clear 
that the diagrams that contribute at order (1/N m , n) are 
related to those that contribute at order (1/N m ~ 1 ,n + 1), 
raising the possibility of a recursive treatment starting 
from the leading order diagrams described here. 

Worth noting is that the semiclassical techniques we 
used here are only valid up to the Heisenberg time, mean- 
ing that we have no access to the density of states on en- 
ergy scales of the order of the mean level spacing. Though 
for ballistic transport the Heisenberg time is much longer 
than the average dwell time (so the mean level spacing 
is much smaller than the Thouless energy) importantly 
the RMT treatment shows that a microscopic gap per- 
sists in this regime even when the time reversal symme- 
try is completely broken (by the magnetic field say). It 
may be possible that applying the semiclassical treat- 
ment of times longer than the Heisenberg time for closed 
system a 79 i 80 to transport would allow one to access this 
regime as well. 

In the opposite regime however, that of the Ehrenfcst 
time, semiclassics provides a surprisingly simple result 49 



allowing complete access to the crossover from the univer- 
sal RMT regime to the more classical Bohr-Sommerfeld 
regime. The gap shrinks due to the suppression of the for- 
mation of encounters while a new class of diagrams (cor- 
related bands) becomes possible. Interestingly the contri- 
bution from trajectories with encounters agrees exactly 
with the results from effective RMT 13 , so our semiclas- 
sical result provides support for this phenomenological 
approach. Of course effective RMT misses the bands of 
correlated trajectories (c.f. those in Ref. liot) which com- 
bined with the other contribution lead to new effects, 
most notably a second gap in the density of states for 
intermediate Ehrenfest times. 
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Appendix: Further generating functions 

The intermediate generating function J(e, r) for the billiard with a single lead and no magnetic field in Sec. UVBl is 
given by 



1 - 



(1-a) 2 + 6r+(l + a)V 1 + 4 (1 - a) 3 - (8 + 20a 2 - a 4 ) r + 4 (1 + a) 3 r 2 
4(1- a) 3 - (16 - 24a + 44a 2 - 8a 3 - a 4 ) r + 2 (12 + 32a 2 - a 4 ) r 2 



- (16 + 24a + 44a 2 + 8a 3 - a 4 ) r 3 +4(1 + a) 3 r 4 



rl a = 0, 



(A.l) 



where we set a = ie. 

The generating function H (e, <p, 6, r) for the billiard with equal leads at phase difference 4> and magnetic field 6 in 
Sec. IV Bl is given by 

f3 2 - ((1 - a + bf + r 2 - 2r (1 - a + 6) (2/3 2 - l)) H 

- r [(1 - a + 6) (1 - 3a + 76) - 2r (l + 56 + b 2 - (3 + 26) a + a 2 ) (2f3 2 - l) + r 2 (1 - 2a + 26)] H 2 
+ r 2 [-6 (196 + 10) + 2a (96 + 1) - 3a 2 + 2r (26 (36 + 4) - 2a (46 + 1) + 2a 2 ) (2/3 2 - l) 

+ r 2 (-6 (6 + 6) + 2a (6 + 1) - a 2 )] H 3 

- r 3 [b (256 + 4) - 14a6 + a 2 - 2r (6 (136 + 4) - 10a6 + a 2 ) (2/3 2 - l) + r 2 (6 (56 + 4) - 6a6 + a 2 )] H 4 

- Ar 4 b [46 - a - 2r (36 - a) (2f3 2 - l) + r 2 (26 - a)] H 5 

- 4r 5 b 2 [l + r 2 - 2r (2/3 2 - l)] H 6 = 0, (A.2) 



where we also used a — ie. For the billiard with unequal leads and no magnetic field in Sec. IV C[ the generating 



22 



function H(e, <p, y, r) is given by 

PP* (1 - a) 2 + pp*r 2 - (/3 2 + /3* 2 ) (1 - a) r 

- (1 - a) 4 + r Ufi + /3* f (1 - a 3 ) + (3 (p + /3*) 2 + 2/3/3*) a (a - 1)) 
+ r 2 ((3 (/3 + /3*) 2 - 2/3/3* - 2) a (2 - a) + 2 (1 + /3 + /3*) (1 - /3 - /3*)) 



+ 



+ r 3 ((/3 + /3*) 2 - a ((/3 + /3*) 2 + 2/3/3*)) - r 4 ] if 

(1 - a) 3 (5a - 1) + ((/3 + /3*) 2 (l - 7a - 7a 3 + a 4 ) + (3/3 + 4/3*) (4/3 + 3/3*) a 2 ) r 

+ 2 (1 + /3 + P*) (1 - /3 - P*) (1 - 6a - 2a 3 ) r 2 ~ (l5/3 2 + 15/3* 2 - 14 + 28/3/3*) 

+ ((/3 + /3*) 2 (1 - 5a) + (3^ 2 + 3/3* 2 + 7/3,3*) a 2 ) r 3 + (4a - 1) r 4 

• 2 [2 (1 - a) 2 (2 - 5a) + (/3 + ^*) 2 (4a 3 - 15a 2 + 15a - 4) r 
+ 2 (1 + [3 + /3*) (1 - (3 - /3*) (a 3 - 8a 2 + 12a - 4) r 2 
+ (/3 + /3*) 2 (-3a 2 + 9a - 4) r 3 + (4 - 6a) r 4 | H 3 



*\ 22 
1 a r 



+ a 2 r 3 



16a - 10a 2 - 6 + (/3 + /3*) 2 (6 - 13a + 6a 2 ) r + 2 (1 + (3 + /3*) (1 - /3 - /3*) (6 - 10a + 3a 2 ) r 2 



+ (/3 + (3* f (6 - 7a + a 2 ) r 3 + (4a - 6) r 4 



H 4 



+ a 3 r 4 



4 - 5a + 4 (/3 + /3*) 2 (a - 1) r + 2 (1 + p + p*) (1 - f3 - /3*) (3a - 4) r 2 



+ (/3 + /3*) 2 (2a - 4) r 3 + (4 - a) r 4 
+ a 4 r 5 (-1 - r 4 + r (l + r 2 ) (/3 + /3*) 2 + 2r 2 [1 + /3 + /3*) (1 - /3 - /3*) 
likewise with a = ie. ^ 



H b = 0, 



(A.3) 
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